%----------------------------------------------
% Formatting and file path
%----------------------------------------------

clear;
clc;
format bank; % Set formatting to display numbers without e format

%----------------------------------------------
% Read in data
%----------------------------------------------
Folder = cd; 
Folder = fullfile(Folder, '../../Data/structural/');
cd(Folder); 
data = csvread('intermediate/data_for_matlab_nonstups.csv');
hhid=data(:,1);
l_bl=data(:,2);
h_bl=data(:,3);
tot_hrs_bl=data(:,4);
k_bl_post_trans=data(:,5);
hired_l_bl=data(:,6);
w_bl=data(:,7);
H=data(:,8);
case_bl=data(:,9);

N=length(hhid);    % Number of HHs
R=3650;            % Time endowment constraint

a=-1.88e-10;       % Estimated parameters
b=0.0009368;
beta=0.3499541;

uns_threshold=10380;

phi = 0.7;

%----------------------------------------------
% Simplify terms
%----------------------------------------------

fk_bl = a.*k_bl_post_trans.^2 + b.*k_bl_post_trans;

%----------------------------------------------
% Calibrate psi_l and A parameters
%----------------------------------------------

psi_l = zeros(N,1);
A = zeros(N,1);

for i = 1:N
    A(i) = phi*w_bl(i)/(fk_bl(i)*beta*(l_bl(i) + hired_l_bl(i))^(beta - 1));
    psi_l(i) = phi*w_bl(i)/l_bl(i);
end

%----------------------------------------------
% Find median psi_l for all non-stups
%----------------------------------------------

median_psi_l_rich=median(psi_l);

%----------------------------------------------
% Find median psi_l for those with log(k)>5
%----------------------------------------------

logk=log((k_bl_post_trans./1000)+1);
psi_l_richest=zeros(N,1);
for i=1:N
    if logk(i)>5
        psi_l_richest(i)=psi_l(i);
    end
end
psi_l_richest=psi_l_richest(psi_l_richest~=0);
median_psi_l_richest=median(psi_l_richest);

save('intermediate/Calibrated_A_psil_nonstups')
